package ikelib;
//-----------------------------------------------------
// http://sakura.bb-west.ne.jp/spr/damayan/algo/spline.html
//
/* スプライン補間(１次元) */
public class Spline implements Mapping {

	/** 補助データの作成 */
	private double[] x, y, z;
	private int N;

	// コンストラクタ
	// x.length == y.length であること
	// x[0] < x[1] < ... < x[N-1] であること
	public Spline(double[] x, double[] y) {
		this.x = x;
		this.y = y;
		N = x.length;
		z = new double[N];
		maketable();
	}

	private void maketable() {
		double[] h = new double[N];
		double[] d = new double[N];

		// 両端点での y''(x) / 6 ... (自然スプライン)
		z[0] = 0;  z[N - 1] = 0;
		for(int i=0; i<N-1; i++) {
			h[i    ] =  x[i + 1] - x[i];
			d[i + 1] = (y[i + 1] - y[i]) / h[i];
		}
		z[1] = d[2] - d[1] - h[0] * z[0];
		d[1] = 2 * (x[2] - x[0]);
		for(int i=1; i<N-2; i++) {
			double t = h[i] / d[i];
			z[i + 1] = d[i + 2] - d[i + 1] - z[i] * t;
			d[i + 1] = 2 * (x[i + 2] - x[i]) - h[i] * t;
		}
		z[N - 2] -= h[N - 2] * z[N - 1];
		for(int i=N-2; i>0; i--) {
			z[i] = (z[i] - h[i] * z[i + 1]) / d[i];
		}
	}

	public double map(double t) {
		int i, j, k;
		double d, h;

		i = 0;  j = N - 1;
		while (i < j) {
			k = (i + j) / 2;
			if (x[k] < t) i = k + 1;  else j = k;
		}
		if (i > 0) i--;
		h = x[i + 1] - x[i];  d = t - x[i];
		return (((z[i + 1] - z[i]) * d / h + z[i] * 3) * d
		+ ((y[i + 1] - y[i]) / h
		- (z[i] * 2 + z[i + 1]) * h)) * d + y[i];
	}

	// テスト
	public static void main(String[] args) {
		int n = 11;
		double[] x = new double[n];
		double[] y = new double[n];
		for(int i=0; i<n; i++) {
			x[i] = -5 + i;
			y[i] = 1/(1 + x[i]*x[i]);
		}
		Spline s = new Spline(x, y);
		for(int i=0; i<10*n+1; i++) {
			double u = -5.5+i*0.1;
			double v = s.map(u);
			System.out.println(u+","+v);
		}
	}

}
/*
f(-5.5)=0.0293056736081598
f(-5.0)=0.038461538461538464
f(-4.5)=0.04761740331491712
f(-4.0)=0.058823529411764705
f(-3.5)=0.07468172670683233
f(-3.0)=0.1
f(-2.5)=0.1400810292242694
f(-2.0)=0.2
f(-1.5)=0.29734709757256067
f(-1.0)=0.5
f(-0.5)=0.8205305804854879
f(0.0)=1.0
f(0.5)=0.8205305804854879
f(1.0)=0.5
f(1.5)=0.2973470975725607
f(2.0)=0.20000000000000007
f(2.5)=0.1400810292242694
f(3.0)=0.10000000000000002
f(3.5)=0.07468172670683233
f(4.0)=0.058823529411764705
f(4.5)=0.04761740331491712
f(5.0)=0.038461538461538464
f(5.5)=0.029305673608159794
*/
